第９・１０回課題

R Markdown

chapter 6

library(sf)

## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

library(raster)

## Loading required package: sp

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:raster':
## 
##     intersect, select, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(spData)
library(spDataLarge)

london = data.frame(lon = -0.1, lat = 51.5) %>% #緯度経度指定
  st_as_sf(coords = c("lon", "lat"))
st_is_longlat(london)

## [1] NA

#> [1] NA

london_geo = st_set_crs(london, 4326) #crsを付加(?)
st_is_longlat(london_geo)

## [1] TRUE

london_buff_no_crs = st_buffer(london, dist = 1)
london_buff = st_buffer(london_geo, dist = 1)

## Warning in st_buffer.sfc(st_geometry(x), dist, nQuadSegs, endCapStyle =
## endCapStyle, : st_buffer does not correctly buffer longitude/latitude data

## dist is assumed to be in decimal degrees (arc_degrees).

london_proj = data.frame(x = 530000, y = 180000) %>% #m指定
  st_as_sf(coords = 1:2, crs = 27700) #地球の球面を平面のように扱える直交座標系の1つ(今回のはBritish national grid)

st_crs(london_proj)

## Coordinate Reference System:
##   User input: EPSG:27700 
##   wkt:
## PROJCRS["OSGB 1936 / British National Grid",
##     BASEGEOGCRS["OSGB 1936",
##         DATUM["OSGB 1936",
##             ELLIPSOID["Airy 1830",6377563.396,299.3249646,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4277]],
##     CONVERSION["British National Grid",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",49,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-2,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996012717,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",400000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",-100000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["UK - Britain and UKCS 49°46'N to 61°01'N, 7°33'W to 3°33'E"],
##         BBOX[49.75,-9.2,61.14,2.88]],
##     ID["EPSG",27700]]

london_proj_buff = st_buffer(london_proj, 111320) #緯度経度を1度=111kmとして換算し、半径1度の円を書く
#緯度経度で半径を設定すると楕円になる。m換算でやると円になる。

Buffers around London with a geographic (left) and projected (right) CRS. The gray outline represents the UK coastline.

#緯度経度を入れることでutmの座標系コードが返ってくる関数の例
lonlat2UTM = function(lonlat) {
  utm = (floor((lonlat[1] + 180) / 6) %% 60) + 1
  if(lonlat[2] > 0) {
    utm + 32600
  } else{
    utm + 32700
  }
}

#動作例
epsg_utm_auk = lonlat2UTM(c(174.7, -36.9))
epsg_utm_lnd = lonlat2UTM(st_coordinates(london))
st_crs(epsg_utm_auk)$proj4string

## [1] "+proj=utm +zone=60 +south +datum=WGS84 +units=m +no_defs"

#> [1] "+proj=utm +zone=60 +south +datum=WGS84 +units=m +no_defs"
st_crs(epsg_utm_lnd)$proj4string

## [1] "+proj=utm +zone=30 +datum=WGS84 +units=m +no_defs"

#> [1] "+proj=utm +zone=30 +datum=WGS84 +units=m +no_defs"

crs_lnd = st_crs(cycle_hire_osm)
class(crs_lnd)

## [1] "crs"

#epsg番号を持ってこれる
crs_lnd$epsg

## [1] 4326

#モルワイデ図法
world_mollweide = st_transform(world, crs = "+proj=moll")

Mollweide projection of the world.

world_laea1 = st_transform(world, 
                           crs = "+proj=laea +x_0=0 +y_0=0 +lon_0=0 +lat_0=0")

Lambert azimuthal equal-area projection centered on longitude and latitude of 0.

cat_raster = raster(system.file("raster/nlcd2011.tif", package = "spDataLarge"))

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on GRS80 ellipsoid in Proj4
## definition

crs(cat_raster)

## Coordinate Reference System:
## Deprecated Proj.4 representation:
##  +proj=utm +zone=12 +ellps=GRS80 +units=m +no_defs 
## WKT2 2019 representation:
## PROJCRS["unknown",
##     BASEGEOGCRS["unknown",
##         DATUM["Unknown based on GRS80 ellipsoid",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1],
##                 ID["EPSG",7019]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8901]]],
##     CONVERSION["UTM zone 12N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-111,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]],
##         ID["EPSG",16012]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]]]

unique(cat_raster)

##  [1] 11 21 22 23 31 41 42 43 52 71 81 82 90 95

wgs84 = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
cat_raster_wgs84 = projectRaster(cat_raster, crs = wgs84, method = "ngb") #この関数でラスター系に変換

plot(cat_raster)

plot(cat_raster_wgs84) #上手く変換できずNULL値の場所があるのが分かる

Key attributes in the original (‘cat_raster’) and projected (‘cat_raster_wgs84’) categorical raster datasets.

CRS	nrow	ncol	ncell	resolution	unique_categories
NAD83	1359	1073	1458207	31.5275	14
WGS84	1394	1111	1548734	0.0003	15

chapter 7

#外部からのデータ取り込み
library(sf)
library(raster)
library(dplyr)
library(spData)

#以下はリンク切れしている。このようにデータをダウンロードする。
#download.file(url = "http://nrdata.nps.gov/programs/lands/nps_boundary.zip",
#              destfile = "nps_boundary.zip")
#unzip(zipfile = "nps_boundary.zip")
#usa_parks = st_read(dsn = "nps_boundary.shp")
#ベクタファイルはst_read

library(rnaturalearth)
usa = ne_countries(country = "United States of America") # United States borders
class(usa)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

#これは古いクラス。こういったものはsfクラスに変換して扱うのが良い。
usa_sf = st_as_sf(usa)

#オープンストリートマップ
library(osmdata)

## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright

parks = opq(bbox = "leeds uk") %>% 
  add_osm_feature(key = "leisure", value = "park") %>% 
  osmdata_sf()

#chapter7.5 file formats

Selected spatial file formats.

Name	Extension	Info	Type	Model
ESRI Shapefile	.shp (the main file)	Popular format consisting of at least three files. No support for: files > 2GB; mixed types; names > 10 chars; cols > 255.	Vector	Partially open
GeoJSON	.geojson	Extends the JSON exchange format by including a subset of the simple feature representation.	Vector	Open
KML	.kml	XML-based format for spatial visualization, developed for use with Google Earth. Zipped KML file forms the KMZ format.	Vector	Open
GPX	.gpx	XML schema created for exchange of GPS data.	Vector	Open
GeoTIFF	.tif/.tiff	Popular raster format. A TIFF file containing additional spatial metadata.	Raster	Open
Arc ASCII	.asc	Text format where the first six lines represent the raster header, followed by the raster cell values arranged in rows and columns.	Raster	Open
R-raster	.gri, .grd	Native raster format of the R-package raster.	Raster	Open
SQLite/SpatiaLite	.sqlite	Standalone relational database, SpatiaLite is the spatial extension of SQLite.	Vector and raster	Open
ESRI FileGDB	.gdb	Spatial and nonspatial objects created by ArcGIS. Allows: multiple feature classes; topology. Limited support from GDAL.	Vector and raster	Proprietary
GeoPackage	.gpkg	Lightweight database container based on SQLite allowing an easy and platform-independent exchange of geodata	Vector and raster	Open

sf_drivers = st_drivers()　#パッケージ読み込みドライバを色々読み込む
head(sf_drivers, n = 5)

##                    name                                  long_name write  copy
## PCIDSK           PCIDSK                       PCIDSK Database File  TRUE FALSE
## netCDF           netCDF                 Network Common Data Format  TRUE  TRUE
## PDS4               PDS4               NASA Planetary Data System 4  TRUE  TRUE
## JP2OpenJPEG JP2OpenJPEG JPEG-2000 driver based on OpenJPEG library FALSE  TRUE
## PDF                 PDF                             Geospatial PDF  TRUE  TRUE
##             is_raster is_vector  vsi
## PCIDSK           TRUE      TRUE TRUE
## netCDF           TRUE      TRUE TRUE
## PDS4             TRUE      TRUE TRUE
## JP2OpenJPEG      TRUE      TRUE TRUE
## PDF              TRUE      TRUE TRUE

vector_filepath = system.file("shapes/world.gpkg", package = "spData")
world = st_read(vector_filepath)

## Reading layer `world' from data source `/usr/local/lib/R/site-library/spData/shapes/world.gpkg' using driver `GPKG'
## Simple feature collection with 177 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## Geodetic CRS:  WGS 84

cycle_hire_txt = system.file("misc/cycle_hire_xy.csv", package = "spData")
cycle_hire_xy = st_read(cycle_hire_txt, options = c("X_POSSIBLE_NAMES=X",
                                                    "Y_POSSIBLE_NAMES=Y"))

## options:        X_POSSIBLE_NAMES=X Y_POSSIBLE_NAMES=Y 
## Reading layer `cycle_hire_xy' from data source `/usr/local/lib/R/site-library/spData/misc/cycle_hire_xy.csv' using driver `CSV'
## Simple feature collection with 742 features and 7 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -0.2367699 ymin: 51.45475 xmax: -0.002275 ymax: 51.54214
## CRS:           NA

#well known text(well known binary)はデータ圧縮率を高めるための形式
world_txt = system.file("misc/world_wkt.csv", package = "spData")
world_wkt = read_sf(world_txt, options = "GEOM_POSSIBLE_NAMES=WKT")
# the same as
world_wkt = st_read(world_txt, options = "GEOM_POSSIBLE_NAMES=WKT", 
                    quiet = TRUE, stringsAsFactors = FALSE, as_tibble = TRUE)

#ラスターデータ
raster_filepath = system.file("raster/srtm.tif", package = "spDataLarge")
single_layer = raster(raster_filepath)

multilayer_filepath = system.file("raster/landsat.tif", package = "spDataLarge")
band3 = raster(multilayer_filepath, band = 3)

multilayer_brick = brick(multilayer_filepath) #メモリに乗せる
multilayer_stack = stack(multilayer_filepath) #メモリに乗せず、ポインタのみ
plot(multilayer_brick)

# data output
st_write(obj = world, dsn = "world.gpkg")

## Writing layer `world' to data source `world.gpkg' using driver `GPKG'
## Writing 177 features with 10 fields and geometry type Multi Polygon.

st_write(obj = world, dsn = "world.shp") #.shpは面倒臭い(大量のメッセージファイル、複数のファイルで1つのデータ)

## Writing layer `world' to data source `world.shp' using driver `ESRI Shapefile'
## Writing 177 features with 10 fields and geometry type Multi Polygon.

## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: Value 318622525 of field pop
## of feature 4 not successfully written. Possibly due to too larger number with
## respect to field width

## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: Value 255131116 of field pop
## of feature 8 not successfully written. Possibly due to too larger number with
## respect to field width

## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: Value 143819666 of field pop
## of feature 18 not successfully written. Possibly due to too larger number with
## respect to field width

## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: Value 124221600 of field pop
## of feature 27 not successfully written. Possibly due to too larger number with
## respect to field width

## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: Value 204213133 of field pop
## of feature 29 not successfully written. Possibly due to too larger number with
## respect to field width

## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: Value 176460502 of field pop
## of feature 56 not successfully written. Possibly due to too larger number with
## respect to field width

## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: Value 1293859294 of field pop
## of feature 98 not successfully written. Possibly due to too larger number with
## respect to field width

## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: Value 159405279 of field pop
## of feature 99 not successfully written. Possibly due to too larger number with
## respect to field width

## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: Value 185546257 of field pop
## of feature 102 not successfully written. Possibly due to too larger number with
## respect to field width

## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: Value 1364270000 of field pop
## of feature 139 not successfully written. Possibly due to too larger number with
## respect to field width

## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: Value 100102249 of field pop
## of feature 147 not successfully written. Possibly due to too larger number with
## respect to field width

## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: Value 127276000 of field pop
## of feature 155 not successfully written. Possibly due to too larger number with
## respect to field width

st_write(obj = world, dsn = "world.gpkg", append = TRUE) #既存ファイルへデータ追加

## Updating layer `world' to data source `world.gpkg' using driver `GPKG'
## Updating existing layer world
## Writing 177 features with 10 fields and geometry type Multi Polygon.

Data types supported by the raster package.

Data type	Minimum value	Maximum value
LOG1S	FALSE (0)	TRUE (1)
INT1S	-127	127
INT1U	0	255
INT2S	-32,767	32,767
INT2U	0	65,534
INT4S	-2,147,483,647	2,147,483,647
INT4U	0	4,294,967,296
FLT4S	-3.4e+38	3.4e+38
FLT8S	-1.7e+308	1.7e+308

writeRaster(single_layer, filename = "my_raster.tif", datatype = "INT2U") #ラスター版出力
writeRaster(x = single_layer,
            filename = "my_raster2.tif",
            datatype = "INT2U",
            options = c("COMPRESS=DEFLATE"),
            overwrite = TRUE) #細かく指定

#image output
#png(filename = "lifeExp.png", width = 500, height = 350)
#plot(world["lifeExp"])
#dev.off()

#こっちでも可
#library(tmap)
#tmap_obj = tm_shape(world) + tm_polygons(col = "lifeExp")
#tmap_save(tm = tmap_obj, filename = "lifeExp_tmap.png")

# htmlなら動く図なんかも保存可能(らしい)
#library(mapview)
#mapview_obj = mapview(world, zcol = "lifeExp", legend = TRUE)
#mapshot(mapview_obj, file = "my_interactive_map.html")

#宿題
#administrative area
unzip(zipfile = "../data/JPN_adm.zip", exdir="./JPN_adm")
jpn_adm0 = st_read(dsn = "./JPN_adm/JPN_adm0.shp")

## Reading layer `JPN_adm0' from data source `/home/rstudio/class/l09-10/JPN_adm/JPN_adm0.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 70 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 122.9332 ymin: 24.04542 xmax: 153.9869 ymax: 45.52279
## Geodetic CRS:  WGS 84

plot(jpn_adm0)

## Warning: plotting the first 9 out of 70 attributes; use max.plot = 70 to plot
## all

#elevation
unzip(zipfile = "../data/JPN_alt.zip", exdir="./JPN_alt")
jpn_alt = raster("./JPN_alt/JPN_alt.gri")

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on WGS84 ellipsoid in Proj4
## definition

plot(jpn_alt)

#population
unzip(zipfile = "../data/JPN_pop.zip", exdir="./JPN_pop")
jpn_pop = raster("./JPN_pop/JPN_pop.gri")
plot(jpn_pop)

jpn_adm0

## Simple feature collection with 1 feature and 70 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 122.9332 ymin: 24.04542 xmax: 153.9869 ymax: 45.52279
## Geodetic CRS:  WGS 84
##   ID_0 ISO NAME_0 OBJECTID_1 ISO3 NAME_ENGLI NAME_ISO NAME_FAO NAME_LOCAL
## 1  114 JPN  Japan        118  JPN      Japan    JAPAN    Japan     Nippon
##   NAME_OBSOL           NAME_VARIA NAME_NONLA NAME_FRENC NAME_SPANI NAME_RUSSI
## 1       <NA> Japan|Ryukyu Islands       <NA>      Japon      Japón     Япония
##   NAME_ARABI NAME_CHINE WASPARTOF CONTAINS SOVEREIGN ISO2  WWW FIPS ISON
## 1    اليابان       日本      <NA>     <NA>     Japan   JP <NA>   JA  392
##    VALIDFR VALIDTO   POP2000     SQKM  POPSQKM    UNREGION1 UNREGION2
## 1 19470402 Present 127096314 373049.4 340.6957 Eastern Asia      Asia
##   DEVELOPING CIS Transition OECD WBREGION          WBINCOME              WBDEBT
## 1          2   0          0    1     <NA> High income: OECD Debt not classified
##   WBOTHER CEEAC CEMAC CEPLG COMESA EAC ECOWAS IGAD IOC MRU SACU UEMOA UMA PALOP
## 1    <NA>     0     0     0      0   0      0    0   0   0    0     0   0     0
##   PARTA CACM EurAsEC Agadir SAARC ASEAN NAFTA GCC CSN CARICOM EU CAN ACP
## 1     0    0       0      0     0     0     0   0   0       0  0   0   0
##   Landlocked AOSIS SIDS Islands LDC                       geometry
## 1          0     0    0       1   0 MULTIPOLYGON (((123.7879 24...

chapter 8

library(sf)
library(raster)
library(dplyr)
library(spData)
library(spDataLarge)

library(tmap)    # for static and interactive maps
library(leaflet) # for interactive maps
library(ggplot2) # tidyverse data visualization package

#tmap(描画ライブラリ)の基本的な使い方

# Add fill layer to nz shape
tm_shape(nz) +
  tm_fill() 

# Add border layer to nz shape
tm_shape(nz) +
  tm_borders() 

# Add fill and border layers to nz shape
tm_shape(nz) +
  tm_fill() +
  tm_borders()

map_nz = tm_shape(nz) + tm_polygons()
map_nz1 = map_nz +
  tm_shape(nz_elev) + tm_raster(alpha = 0.7) #重ねているが、透明度を指定しているので透けて見える
nz_water = st_union(nz) %>% st_buffer(22200) %>% 
  st_cast(to = "LINESTRING")
map_nz2 = map_nz1 +
  tm_shape(nz_water) + tm_lines()
map_nz3 = map_nz2 +
  tm_shape(nz_height) + tm_dots()
tmap_arrange(map_nz1, map_nz2, map_nz3)

## stars object downsampled to 877 by 1140 cells. See tm_shape manual (argument raster.downsample)
## stars object downsampled to 877 by 1140 cells. See tm_shape manual (argument raster.downsample)
## stars object downsampled to 877 by 1140 cells. See tm_shape manual (argument raster.downsample)
## stars object downsampled to 877 by 1140 cells. See tm_shape manual (argument raster.downsample)

## Some legend labels were too wide. These labels have been resized to 0.62, 0.54, 0.54, 0.54, 0.54, 0.54. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

## stars object downsampled to 877 by 1140 cells. See tm_shape manual (argument raster.downsample)

## Some legend labels were too wide. These labels have been resized to 0.62, 0.54, 0.54, 0.54, 0.54, 0.54. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

## stars object downsampled to 877 by 1140 cells. See tm_shape manual (argument raster.downsample)

## Some legend labels were too wide. These labels have been resized to 0.62, 0.54, 0.54, 0.54, 0.54, 0.54. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

ma1 = tm_shape(nz) + tm_fill(col = "red")
ma2 = tm_shape(nz) + tm_fill(col = "red", alpha = 0.3)
ma3 = tm_shape(nz) + tm_borders(col = "blue")
ma4 = tm_shape(nz) + tm_borders(lwd = 3)
ma5 = tm_shape(nz) + tm_borders(lty = 2)
ma6 = tm_shape(nz) + tm_fill(col = "red", alpha = 0.3) +
  tm_borders(col = "blue", lwd = 3, lty = 2)
tmap_arrange(ma1, ma2, ma3, ma4, ma5, ma6)

#plotとtmapの違い
plot(nz["Land_area"])

tm_shape(nz) + tm_fill(col = "Land_area")

legend_title = expression("Area (km"^2*")") #タイトルの付け方
map_nza = tm_shape(nz) +
  tm_fill(col = "Land_area", title = legend_title) + tm_borders()

tm_shape(nz) + tm_polygons(col = "Median_income") #デフォルトは6段階に分ける

breaks = c(0, 3, 4, 5) * 10000
tm_shape(nz) + tm_polygons(col = "Median_income", breaks = breaks) #0,3万,4万,5万で分ける

tm_shape(nz) + tm_polygons(col = "Median_income", n = 10) #10段階で分ける

tm_shape(nz) + tm_polygons(col = "Median_income", palette = "BuGn") #色を変える

tm_shape(nz) + tm_polygons("Population", palette = "Blues")

## Some legend labels were too wide. These labels have been resized to 0.63, 0.58, 0.58. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

tm_shape(nz) + tm_polygons("Population", palette = "YlOrBr")

## Some legend labels were too wide. These labels have been resized to 0.63, 0.58, 0.58. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

map_nz + 
  tm_compass(type = "8star", position = c("left", "top")) +
  tm_scale_bar(breaks = c(0, 100, 200), text.size = 1)

#色々オシャレにできる
map_nz + tm_layout(title = "New Zealand")

map_nz + tm_layout(scale = 5)

map_nz + tm_layout(bg.color = "lightblue")

map_nz + tm_layout(frame = FALSE)